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We show that spatial models of simple predator-prey interactions predict that predator and prey 
numbers oscillate in time and space. These oscillations are not seen in the deterministic versions 
of the models, but are due to stochastic fluctuations about the time-independent solutions of the 
deterministic equations which are amplified due to the existence of a resonance. We calculate the 
power spectra of the fluctuations analytically and show that they agree well with results obtained 
from stochastic simulations. This work extends the analysis of these quasi-cycles from that previ- 
ously developed for well-mixed systems to spatial systems, and shows that the ideas and methods 
used for non-spatial models naturally generalize to the spatial case. 



PACS numbers: 87.23.Cc,02.50.Ey,05.40.-a 
I. INTRODUCTION 

The standard paradigm of condensed matter physics 
involves the interaction of discrete entities (for example 
atoms, molecules or spins) positioned on the sites of a reg- 
ular lattice which when viewed at the macroscale can be 
described by a differential equation after coarse-graining. 
This type of structure is not unique to physics; there 
are many other systems which consist of a large num- 
ber of discrete entities which interact with each other in 
a simple way, but which when viewed macroscopically 
show complex behavior. What is different, however, is 
that physicists stress the relationships between models 
of the same phenomena constructed at different scales, 
for instance by deriving macroscopic models from those 
defined at the microscale. Here we will be interested in 
modeling species in an ecological system where the in- 
teraction between individuals of those species is of the 
predator-prey type. Although both "microscopic mod- 
els" — individual based models (IBMs) defined on a 
two-dimensional lattice for example, and "macroscopic 
models" such as reaction-diffusion equations, have been 
extensively studied , the derivation of the latter from 
the former has received very little attention. Thus it is 
not obvious a priori if the results from the two differ- 
ent approaches can be meaningfully compared or if the 
macroscopic description misses some important features 
which are present in the IBM. 

In this paper we will build on some earlier work 0] that 
introduced a methodology which began from a specific 
IBM and derived the corresponding model which holds 
at the macroscopic, or population, level. The latter was 
called the population level model (PLM) and the for- 
mer sometimes called the individual level model (ILM), 
rather than the IBM, by analogy. There is another reason 
for using the term ILM in place of IBM. The nature of 
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the "microscopic model" can vary considerably. At one 
extreme are models where the constituents each have in- 
dividual characteristics. They may have an age, sex, be 
hungry at a given time, and so on. These are essentially 
agent based models [!, [j] . At the other extreme are very 
simple "physical" models, such as lattice gas models [H], 
where the analogies to physical processes take a primary 
role. The term IBM is frequently used for the former 
agent based models. Our starting point will be some- 
where between these two extremes. We model the indi- 
viduals as entities which may be born or die, may migrate 
to neighboring sites on the lattice in a single time interval 
and when on the same lattice site may interact with each 
other if one is a prey species of the other. Thus the in- 
dividuals act as chemical species which have given inter- 
action rules. There are several advantages with this for- 
mulation. Firstly, it corresponds most directly in terms 
of properties of the constituents to PLMs such as the 
Volterra equations. Secondly, more properties of indi- 
viduals can be included if required, taking the model 
more towards the agent-based IBMs mentioned above. 
Thirdly, it allows the stochastic nature of birth, death, 
predator-prey and migrationary processes to be naturally 
included into the model. Whereas most stochastic mod- 
els have been simulated directly, we prefer to formulate 
them as a master equation, and use the system-size ex- 
pansion Q to derive the form they take when the system 
size is large. 

The aim of this paper is then to investigate the nature 
of the PLM model both at the macroscopic or mean-field 
level — which is deterministic — and at what might be 
described as the mesoscopic scale where stochastic effects 
are still important, but where the discrete nature of the 
constituents has been lost. The former is interesting be- 
cause it is not clear that the model derived in this way 
will coincide with those appearing in the textbooks on 
the subject 0, H, Ho| . but also because of the types 
of collective patterns frequently displayed by these sys- 
tems, which often resemble those observed when study- 
ing physical and chemical systems. The latter is inter- 
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esting because it has been found that in simple predator- 
prey models (without spatial effects being included) large 
predator-prey cycles are present in the stochastic model, 
which are lost at the deterministic level [lf[ . More specifi- 
cally, the discrete nature of the individuals results in a de- 
mographic stochasticity at the mesoscale which acts as a 
driving force and creates a resonance effect, turning small 
cyclic fluctuations into large cycles called quasi-cycles [8| . 
Here we investigate the nature of this phenomenon in a 
model where spatial effects are included. The ordinary 
differential equations of the Volterra type will now be 
replaced by partial differential equations of the reaction- 
diffusion type, and the two coupled Langevin equations 
of [ll| will be replaced by two coupled partial differential 
equations with additive noise. 

The paper is organized as follows. In Section HH the 
model alluded to above is introduced and formulated as a 
master equation. This is followed in Section 111 by a dis- 
cussion of the deterministic limit of the equation, a linear 
stability analysis of the stationary solution of this equa- 
tion, and the linear noise correction to the deterministic 
equation. In Section IV a Fourier analysis of the lin- 
ear stochastic differential equations is carried out which 
yields power spectra which characterize the nature of the 
spatial and temporal predator-prey cycles, with the an- 
alytic results being compared to the results of computer 
simulations. There are two Appendices containing math- 
ematical details. The first describes the application of the 
system-size expansion to the master equation and the sec- 
ond contains the Fourier analysis of the linear stochastic 
differential equation. 



II. MODEL 



are assumed to be local, that is, only involve individuals 
at a particular site. They will therefore be identical to 
those invoked in the predator-prey model without spatial 
structure [ll[ , and since these have been shown to lead 
to the Volterra equations in the deterministic limit, we 
will adopt them here: 
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All constituents have a subscript i to denote that they 
are located in patch i. Eq. ([I]) describes the birth of a 
prey individual, which occurs at a rate b. We assume 
that "space" is required for this to occur. Also we do 
not specify the birth of predator individuals as a sepa- 
rate event, since these also occur through predation, as 
described by Eq. and will not lead to new terms in 
the evolution equations. Two types of predation are re- 
quired in Eq. ||2J| so that only a fraction of the resources 
obtained from consumption of the prey are used to pro- 
duce new predator individuals. Finally, Eq. Q describes 
the death of individuals of species A and B at rates d\ 
and g?2 respectively. 

Here we arc considering an explicitly spatial model, so 
the additional feature which we include is the possibility 
of changes in the populations due to migrations between 
nearest neighbor patches. These events can be described 
by adding the following set of reactions Q : 



The system we will be interested in consists of individ- 
uals of species A who are predators of individuals belong- 
ing to the prey species B. We assume that they inhabit 
patches, labeled by i = 1, . . . , f2, which are situated at the 
sites of a d— dimensional hypercubic lattice. Of course, 
for applications we are interested in the case of a square 
lattice in two dimensions, but we prefer to work with gen- 
eral d. One reason is that it is not any more complicated 
doing so, another is because our stochastic simulations 
have been carried out in d = 1 in order to achieve higher 
accuracy. Each patch possesses a finite carrying capacity, 
N, which is the maximum number of individuals allowed 
per site. The number of predators and prey in patch i 
will be denoted by rij and nij respectively. There will 
therefore be (N — — mi) empty or vacant "spaces" , E, 
in patch i. These are necessary to allow the number of A 
and B individuals in patch i to independently vary with 
time. Further background to the modeling procedure is 
given in Ref. 0] , where it has been applied to competition 
between two species. 

As discussed in Section HI we assume that the con- 
stituents A, B and E react together at given rates. The 
reactions corresponding to birth, death and predation 



A%Ej — EiAj ; BiEj EiBj , 

A 3 Ei ^ ; BjEi ^ EjBi . (4) 

Here i and j are nearest neighbor sites and /ii and \ii 
are the migration rates for individuals of species A and 
B respectively. 

The state of the system at any given time is speci- 
fied by the elements of the set {n.;, m, : i = 1, . . . , f]}. If 
we take the transition rates between these states to only 
depend on the current state of the system, the process 
will be Markov and can be described by a master equa- 
tion in continuous time. The natural way to define such 
transition rates is according to a mass action law: the 
probability that two constituents meet is proportional to 
their current proportions in their respective patches. The 
allowed transitions and the rates at which they take place 
are given by Eqs. (HJ-Jl]). Denoting the transition rates 
from a state with ni predators and m,k prey to a state 
with n\ predators and m' k prey by T n / t7n > \ numk , then the 
transition rates corresponding to the purely local reac- 
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tions d)-© are: 

-^ri.i+l,mj— l\ni,mi 



function of n and m as follows: 
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These are exactly as in the non-spatial form of the 
model [UJ , but with the state variables all having a sub- 
script i to denote these are the reactions in patch i and 
an extra factor of Q in the denominator since there is 
a choice between any one of the O patches when deter- 
mining the probability of a transition taking place. To 
lighten the notation we have shown the dependence of T 
only on the subset of variables liable to change (in this 
case those on the site i). The corresponding expressions 
for the transition rates between nearest neighbors, which 
describes the migratory process, are 
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Here, z denotes the coordination number of the lattice, 
that is the number of nearest neighbors of any given site, 
which in our case is 2d. It needs to be included since 
it represents the choice of nearest neighbor j, once the 
patch i has been chosen. 

The master equation which governs the time evolu- 
tion of the system can now be constructed. Although 
this equation can easily be written down, and has the 
standard form of a sum of transition probabilities giv- 
ing rise to a change in the probability distribution func- 
tion with time [6[ , it has a rather ungainly appear- 
ance. It can be made to look neater through the in- 
troduction of a little more notation. First, the prob- 
ability distribution function that the system is in state 
{rii, rrii : i = 1, . . . , f2} at time t is conventionally denoted 
by P (m, mi, . . . ,riQ, ma; t), but we will denote it by 
Pn,m(t)- Then the master equation takes the form 



dP n , m (t) 
dt 
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(7) 

where the notation j € i means that j is a nearest neigh- 
bor of i and where 7" loc and T,f ls are transition rates 
which are defined below. These transition rates may in 
turn be simplified by the introduction of the step opera- 
tors @ E^ 1 and E^ 1 defined by their effect on a typical 
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The local transition operator 7{ oc may now be written 
as 
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with the four local transition rates given explicitly in 
Eq. Similarly, the transition operator 7^™ Ig which 

involves transitions between nearest neighbor sites can 
be written as 
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The master equation ([7]), together with the definitions 
of the transitions rates given by Eqs. (© and (O together 
with Eqs. (© and (jTUJ) , completely define the model once 
initial and boundary conditions are specified. The model 
is far too complicated to be solved exactly, but it can be 
analyzed very accurately by studying it in the limit of 
large system size. As previously proposed @, [TJ, H3.ll3j|. 
and as discussed in Appendix A, the leading order in a 
system-size expansion of the master equation gives de- 
terministic equations whose stationary state can be ana- 
lyzed, whereas the next-to-leading order result gives lin- 
ear stochastic differential equations, which can be Fourier 
analyzed. From this we can investigate the possible exis- 
tence of resonant behavior induced by the demographic 
stochasticity of the original model. 

In the next section we analyze the equations describing 
the model to leading order and next-to-leading order in 
the system-size expansion. The details of the calculation 
required to determine these is given in Appendix A. 



III. DETERMINISTIC LIMIT AND 
FLUCTUATIONS ABOUT IT 

The deterministic limit of the model defined by 
Eqs. 0, © and Jl0|) is derived in Appendix A. It is 
defined in terms of the populations <pi — Ym\N^oo{ni/N) 
and ipi = limAr^ 00 (mi/iV) and explicitly given by Eqs. 
EH), (lA"5j) . (|A15|) and (IA16|) . These may be written as 
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the 2Q macroscopic equations 
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where i = 1, . . . , Q and where the symbol A represents 
the discrete Laplacian operator Afc — - (fj ~ /*)■ 
A rescaled time, r = t/Q, has also been introduced. 

To complete the formulation of the problem, initial 
and boundary data should be provided. For the type 
of system considered here the most natural choice is to 
consider zero-flux boundary conditions, regardless of the 
initial conditions. This corresponds to the condition that 
individuals are not allowed to leave or enter the fixed re- 
gion designated as the system, in other words there is 
no immigration or emigration. The system of equations 
(fTTj) - (|12[) possesses two limits of interest. The limit = 1 
formally corresponds to a one-site system and is simply 
the well-known Volterra model as studied in [ill ]. The 
limit f2 — ► 00 corresponds to shrinking the lattice spac- 
ing to zero and so obtaining a continuum description in 
which the discrete Laplacian operator is replaced by the 
continuous Laplacian V 2 and the Eqs. (| 1 1 [) - (fl~2|) become 
a pair of partial differential equations: 
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where a = 2p l7 = d 1: r = 2b - d 2 , K = (26 - d 2 )/2b, 
and A = 2 (pi +P2 + b), with ip and <f> representing the 
prey and predators densities respectively. It should be 
noted that in the transition to a continuum model, the 
population fractions go over to population densities and 
parameters may be scaled by factors involving the lat- 
tice spacing. An example of this involves the migration 
rates in Eqs. (fT3"]) - (|14[) . which are scaled versions of those 
appearing in Eqs. (|TI ]) -([T2" ]) (see Eqs. (|B10j) ). 

One of the most interesting features of Eqs. (]13p - (fl"4]) 
is the emergence of cross-diffusive terms of the type 
(■0V 2 — (f)V 2 ip). These types of contributions do not 
usually appear in the heuristically prop osed spatially ex- 
tended predator-prey models [lj, [l5|. However, they 
seem to appear naturally in these types of lattice models, 
and cross-diffusive terms similar to those found here have 
been obtained as the mean-field limit of a set of models 
proposed by Satulovsky [j^| . An inspection of Eqs. (fT3]l - 
(fT4]> leads to the conclusion that they do not reduce to 
a simple reaction-diffusion scheme for any choice of pa- 
rameters, however if zero-flux boundary conditions are 



chosen, this implies that, after a single integration over 
the spatial domain, the contribution of the cross-diffusive 
terms for the solution vanishes: 
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with a similar equation with (j> and ip interchanged. The 
condition (|15[) also occurs if we impose the requirement 
that 4>{ r >t) an d <Xr, t) vanish as r — > 00, instead of the 
zero-flux boundary conditions, which are those typically 
chosen in textbooks [l0( . 

Before discussing the equations which describe the 
stochastic behavior of the system, we will analyze the 
nature of the stationary solutions in the deterministic 
limit. We will be particularly interested in investigating 
the possibility that "diffusion-driven" instabilities may 
occur for the model defined by Eqs. (fTTj) - (fl"2"]) or equiva- 
lently for Eqs. (pjl - pT]) . 



A. Stationary state in the deterministic limit. 

One of the simplest questions one can ask about 
Eqs. HU-fH} or Eqs. p^] ) -([T3 j) concerns the nature of 
the stationary state. It is simple to verify that there 
are two unstable fixed points (describing the null state, 
(ff = t/j* = 0, and a state with no predators, <f>* = 0, 
tp* = i4T)j_and a single coexistence fixed point given by 
(see also [13, HI, HH , for instance) 



aK J 
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(16) 



Finding non-homogeneous stationary state solutions 
would require solving a pair of coupled non-linear dif- 
ferential equations, but we can look for solutions if the 
homogeneous solutions p^]) are unstable to spatially in- 
homogeneous small perturbations. That is, we look for 
solutions of Eqs. (fTT|) - (fT2]) which have the form 



tpi = ip* + ■ 



(17) 



where Uj and Vj are the small perturbations. An exactly 
similar analysis could be carried out on the continuum 
versions (|13[) - (I14I) . but now u and v would be functions of 
r, a vector in the region of interest. Substituting Eq. (fl~7]) 
into Eqs. (fTTj) - (fl"2"]) . and keeping only linear terms in u 
and v gives 
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Here an,ai2,a2i and 022 are the contributions which 
would be found if the perturbation had been assumed 
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to be homogeneous; they are exactly the terms found in 
[lH . namely 



an = atp* — (3; ai% 
fl2i = -AV>*; a 2 2 



' (1-^1 (20) 



We may write Eqs. (fT5|) and Q19p in the unified form 
iij = Auj with Uj = (uj, Vj) T for a given site j. The en- 
tries of the matrix A will be denoted by am, 0^12, ctj^i 
and ai t 22- The solution to Uj = Aiij has the form 



Uj(r) ~ expjz^r + iak.j} . 



(21) 



where a is the lattice spacing and where we have explic- 
itly indicated the vector nature of j and k. The v and k 
must satisfy 



v — an —oti2 

— tt21 V — «22 
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where 



Ok, 11 = an + hi (1 - ip*) A k , 

«k,12 = «12 + Ml^*Ak , 

£*k,21 = 2 1 + jti 2 ^*A k , 

«k,22 = a 22 + M2 (1 - 0*) A k , 



(23) 



and where the discrete Laplacian, A k for a 
d— dimensional hypercubic lattice is (see Appendix 
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The idea that patterns can form due to a diffusion- 
induced instability was first put forward by Turing in 
1952 in connection with his investigation into the ori- 
gins of morphogenesis [20( . More generally, such patterns 
can arise in reaction-diffusion equations where a homoge- 
neous stationary state is stable to homogeneous perturba- 
tions, but where irregularities or stochastic fluctuations 
in real systems can induce local deviations from the spa- 
tially uniform state, which can in turn grow if this state 
is unstable to inhomogeneous perturbations. Since Tur- 
ing's seminal work, the phenomenon has been studied in 
many types of reaction-diffusion system, including spa- 
tial predator-prey models [3, [U, [2^, HI] • In contrast to 
these previous studies, where the reaction-diffusion equa- 
tions we postulated phenomenologically, we have derived 
our equations from a ILM. Moreover they differ from the 
models considered previously because of the existence of 
non-linear diffusive terms. Therefore it is of interest to 
study if the model we have derived allows for the exis- 
tence of Turing patterns. 

We first need to check that the homogeneous station- 
ary state is stable to homogeneous perturbations. A ho- 
mogeneous perturbation means that the Uj and Vj in 



Eq. (fl~7|) are independent of j. This in turn means that 
the terms involving fii and (12 are absent from Eqs. (|18p 
and (TT9")) . Therefore the stability to homogeneous per- 
turbations may be found from Eq. (|22p with the a re- 
placed by the a. Stability is assured if an + 022 < and 
ano 2 2 — 012021 > 0, since these conditions are equivalent 
to asking that the v which are solutions of Eq. (|2"2"|) have 
negative real parts. It is straightforward to check from 
the explicit forms (p~6|) and ([20]) that an — 0, 012 > and 
021, «22 < 0, and so that this is the case. As an aside we 
can also check that for the null state (cf>* = ip* = 0) and 
the state without predators ((f)* = 0, ip* = K), under the 
condition that the fixed point (|16p exists, that 011022 < 
and 012 = 0. Therefore the determinant of the stability 
matrix is negative, and so the eigenvalues are real with 
different signs, and both these states are unstable. 

To get a diffusive instability, we need to investigate 
the solutions (fTT)) which now include the spatial contri- 
butions. For an instability to occur, one of the conditions 
trA k < or det A k > must be violated. From Eq. (J23]) it 
is clear that A k < and so from Eq. (|2"3")l that an < an 
and «22 Q-22 and so that trA k < 0. So the only possi- 
bility for a Turing pattern to arise is if detA k < 0. By 
direct calculation 

detA k = -oi 2 a 2 i - Mi [o2i0* - a 2 2 (1 -ip*)] A k 

- /i 2 ai 2 V'*A k + (ii(i 2 (1 - 0* - ip*) A£ . (25) 

Now all the terms on the right-hand side of Eq. ([23)1 are 
manifestly positive, except the second. However, since 



- <<ii n - >- > = /-<•' ( - 1 



(26) 



and K = 1 — (<i 2 /2&) < 1, then this term is also positive. 
Therefore detA k > and so the homogeneous station- 
ary state is stable to both small homogeneous and small 
inhomogeneous perturbations. 

It has been known for some time that the simple 
reaction-diffusion equations for a predator-prey model 
(i.e. those containing only containing simple diffusive 
terms such as V 2 (f> and V 2- 0) do not lead to diffusive 
instabilities [l(| • We have shown here that the introduc- 
tion of a particular type of cross-diffusive term, which 
has its origins in the ILM formulation, also contains no 
Turing instability. It should be noted that this also holds 
true in the limit of zero lattice spacing where A k is re- 
placed by — k 2 (up to a constant), which is also always 
negative for k ^ 0. This corresponds to using Eqs. (fT3"|) - 
(fT4")) . rather than Eqs. (fTTj) - (|12p . Since, on average, the 
population fractions do not exhibit any form of spatial 
self-organizing structure, the emergence of such struc- 
tures when observing the full dynamical process should 
be understood as an effect due to fluctuations. So we 
now study the next next-to-leading order contributions 
which describe fluctuations around these mean values, 
with the aim of quantifying possible resonant behavior 
in both space and time. 
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B. Fluctuations 

The next-to-leading order in the system size expansion 
gives a Fokker-Planck equation in the 2f2 variables and 
rji, which describe the deviation of the system from the 
mean fields: 

(27) 

The equation itself is derived in Appendix [AJ it is given 
by Eq. (|A6|) with coefficients defined by Eqs. (|A26|) and 
(IA27|) . These coefficients have been evaluated at the 
fixed-point <jf , ip* of the deterministic equations since, as 
explained earlier, we are interested in studying the effect 
of fluctuations on the system once transient solutions of 
the deterministic equations have died away. Rather than 
work with this Fokker-Planck equation, it is more conve- 
nient to use the Langevin equation which it is equivalent 
to. This has the form @, Hi] 

^=A-(C) + A i (r), (28) 

where 

(X i (r)\ j (T')}=B ij 5(r-r'). (29) 

Here Ci = (6, Vi) and K = (\i, K2) with % being the 
constant matrix defined by Eq. (|A27jl . 

The key point here is that the system-size expansion 
to this order yields a function A(C) which is linear in , 
as can be seen from Eq. (|A26[) . It is this linear nature 
of the Langevin equation which is crucial in the analysis 
that follows. To study possible cyclic behavior we will 
require to calculate the power spectrum of the fluctua- 
tions (|27p. and to do this we need to find an equation 
for their temporal Fourier transforms. The linearity of 
the Langevin equation (|28|) means that this is readily 
achieved. The translational invariance of the solutions of 
the deterministic equations, together with the nature of 
the diffusive terms also make it useful to take the spatial 
Fourier transform of Eq. (j28|) . This is discussed in detail 
in Appendix [Bj writing out the two components of the 
equation explicitly it has the form 

-7^ = Qik.ii £k + ak,i2 ??k + Ai k(r) 
dr 

—j— = a k .2i £k + a k . 22 ?7k + A 2 .k(i~) , (30) 
ar 

where the a k are given by Eq. (f2"3"|) and by Eq. (|2T)|) . The 
noise correlators (f2l?|) are now local in k— space: 

(A k (r)A k , (rO) = B k n a d <5 k+ k',o<5(r - r') , (31) 

where i3 k is derived in the Appendices (see Eq. (|B6[) et 



seq.) and is given by 

B k ,n = a d [(di<f>* + 2piW*) 

- 2 Ml 0*(l-0*-V*)A k ] , 

B K2 2 = a d [2biP* (1 - 4>* - $*) + d 2 r 

+ 2(p 1 +p 2 )V*0*-2 M2 V/(l-^-^)A k ] , 

Bk,i2 = Bk,2i = -2a d Pl( f>*r . (32) 

It should be noted that, since A k < 0, the diagonal ele- 
ments of Bk are all positive, as they should be. 

It is interesting to consider what happens in the con- 
tinuum limit a — > 0. For non-zero a, the wave-numbers 
take on values in the interval (—n/a) < hi < (ir/a), but 
this becomes an infinite interval as a — * 0. The wave- 
numbers are still discrete however, due to the finite vol- 
ume (area in two dimensions) of the system; we keep the 
volume f2a rf fixed in the limit, so that — * 00. In the 
limit Qa d ^ic+ic^o goes over to (27r) d <5(k+k') and A k goes 
over to — fc 2 , as long as the migration rates are suitably 
scaled (see Eq. (jBlOj) 1 ). However from Eq. (|3"2")l . it is clear 
that the 2? k vanish in the limit due to the factor of a d . 
This should not be too surprising: since f2 — ■+ 00, the 
number of degrees of freedom of the system is becoming 
infinitely large, and thus we would expect fluctuations to 
vanish. If all the Z? k are zero, the noises A k (r) vanish, 
and therefore so do £k(r) and ?7 k (r). This effect has been 
seen (see (2(| and the references therein) : oscillatory be- 
havior in these types of models persists as long as the 
number of sites remains finite, however it disappears in 
the so-called thermodynamic limit. However, in practice, 
one has to go over to describing the population sizes as 
population densities, rather than pure numbers, in this 
limit. This will involve further rescalings, and depending 
on the exact definition of the model, these fluctuations 
can survive the continuum limit. For this reason we will 
keep a finite lattice spacing: the results for a particu- 
lar continuum model variant can then be determined by 
taking the a — > limit in the appropriate manner. 



C. Simulations 

We expect that the deterministic equations (TTTj) and 
(TT2")) . together with the stochastic fluctuations about 
them, given by Eqs. (|3"0|) - (f3"2"|) , will give an excellent de- 
scription of the model defined by Eqs. ©-(Ill) for moder- 
ate to large system size. We test this expectation here by 
presenting the results of numerical simulations performed 
for the full stochastic process (H|)-(|ll) using the Gillespie 
algorithm [27J ■ This is completely equivalent to solving 
the full master equation ([7]) . To obtain the best results we 
restricted our simulations to the one-dimensional system 
(d = 1), even though our theoretical treatment applies 
to general d and we would usually be interested in d = 2. 
We took the length of the spatial interval to be unity, so 
that ail = 1. Therefore once the number of lattice sites, 
£7, is fixed, so is the lattice spacing, a. 
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FIG. 1: Results obtained from simulating the process CEJ-©. Panel (a) shows the temporal evolution of the total population 
fractions of predators and prey ^(t). Panels (b)-(f) contain snapshots at different times of the spatial configuration for 
a typical realization (panels (c) and (e)) and averaging 150 independent realizations (panels (d) and (f)). Panel (b) shows 
the initial spatial configuration. The reaction rates employed were pi = 0.25J1, P2 = 0.05S1, di = O.lfi, c?2 = 0.0, b — 0.1S1, 
[ii = 0.2S1, [12 = O.lfi, fi = 200 and N = 500. The dotted lines in the figure correspond to the fixed point values <j>* and ip* 
found from Eq. (fT6)l . 
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Panel (a) of Fig. Q] shows typical behavior of the to- 
tal population fractions <&(i) = Si=i Ui anc ^ = 
TW SiLi m « starting from the initial condition shown in 
Fig. [Ub). Subsequent panels show the time evolution 
of the local fractions <pi an d ipi starting from the same 
initial condition. The time t corresponds to the Gille- 
spie time and was measured in integer time-steps. The 
average values are those calculated from the fixed point 
(fT6]) . For this simulation the number of sites employed 
was 12 = 200 and the site capacity was N = 500. The 
local reaction rates were chosen so as to match the values 
used in the non-spatial version of the model [ll[ . In par- 
ticular this means that <fi* — ip*. Since the time in this 
spatial version is scaled by 12 (t — t/fl), the rates are 12 
times those used in [ll|, namely p\ = 0.2517, p 2 — 0.0512, 
di = 0.1Q, d 2 = 0.0 and b = 0.117. The values of the 
migration rates /ii and p, 2 for this simulation were 0.212 
and O.lfi respectively. 

The initial configuration shown in Fig. QJb) consists of 
prey homogeneously distributed along the spatial inter- 
val, with populations equal to the equilibrium coexistence 
value mi = ip*N. The predator species were also initially 
homogeneously distributed, with the difference that they 
were confined to only the middle third of the sites; the 
first and last third of the interval contained no preda- 
tor individuals. This choice was made in order to clearly 
indicate the nature of the invasion process of predators 
into the predator-free zones, which eventually leads to 
the establishment of a mixed predator-prey regime over 
the whole spatial interval. Before this happens, all those 
sites with only prey individuals should converge to the 
saturation value iji = K and remain there until a preda- 
tor invades the site, which can only occur via a migration 
event. Once the entire spatial interval is populated with 
individuals of the two species, their numbers will oscillate 
around the fixed point (4>*,tp*), as shown in Figs. [TJc) 
and (e) . It was found that for the parameter values taken 
in this realization of the process, the mixed state first be- 
comes established in the entire domain at approximately 
t ~ 800. For times larger than this there is oscillatory 
behavior around the fixed point values, which is shown in 
later figures (Figs.[2[a) and (b)); this behavior resembles 
that reported in [111 ] . 

Figure [1] also contains a sequence which shows the dy- 
namics of the average values (panels (d) and (f)), ob- 
tained by averaging over 150 independent realizations 
of the stochastic process. The dynamics consists of a 
continuous transition from the unstable state with only 
prey present, into the stable two-species fixed point. This 
takes the form of traveling wave-fronts of "pursuit" and 
"evasion" , which describe the invasion process of preda- 
tors into locations occupied only by prey individuals. 
Such traveling waves may be found directly as solutions 
of the deterministic equations [13, EH • 

Our main interest in this paper is the study of the 
nature of the fluctuations about the stationary state, that 
is, at times subsequent to that illustrated in Fig. [ljd), 
and we now return to their study. 



IV. POWER SPECTRA 

To calculate the power spectra of the fluctuations 
about the stationary state, we first have to take the tem- 
poral Fourier transform of Eqs. (|30| . This reduces the 
equations governing the stochastic behavior of the system 
to two coupled algebraic equations which are linear, and 
so which can be used to obtain a closed form expression 
for the power spectra. In this section we first describe 
this analytic approach, and then go on to discuss how 
the power spectra can be found from numerical simula- 
tions, and then finally compared the results of these two 
approaches. 



A. Analytic form 

Taking the temporal Fourier transform of Eqs. (|30|) 
yields 



MC k M = A k H , 



(33) 



where M = (— iul — A) and I is the 2x2 identity matrix. 
Therefore C = M _1 A, which implies that 

|£k(w)| 2 = |pn| 2 A1A1 + p 11 p* 12 X 1 X 2 

+ PI1P12KM + \P2-1\ A2A2 , (34) 

with a similar expression for |?7k(^)| 2 which is just 
Eq. (|3"4"|) but with all the first subscripts of p changed 
to 2. Here the p a b are the elements of M" 1 . Using 



(A k (w)A£H) =B k , 
the power spectra for the predators 



a,iM = (i&Mi 



and for the prey 



k,2 



(35) 



(36) 



(37) 



may easily be found. 

Since the Langevin equations are diagonal in k— space, 
the structure of the expressions for the power spectra 
are the same as those found in other studies [lj], [12, [3 1 
namely 



flu 



C k .i +Buhuj 2 



(^ 2 -^o)' 



(38) 



and 



fl.,2 = 



Ck,2 + $k,22W 2 



(39) 
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where 

Ck,l = #k,ll "k,22 — 2Sk,12 ak,12 «k,22 + #k,22 «k : 12 ) 
Ck,2 = #k,22 a k,ll — 2Sk,12 ak,21 ak.ll + #k,ll «k;,21 ■ 

'(40) 

The spectra ([55)1 and ([3"9"]) resemble those found when 
analyzing driven damped linear oscillators in physical 
systems. A difference between that situation and the 
one here is that the driving forces here are white noises 
A(t) which excite all frequencies equally, thus there is 
no need to tune the frequency of the "driving force" to 
achieve resonance. The parameters in the denominators 
of Eqs. ([38]) and ([39]) are given by fi 2 . = detAk and 
Tk = — trAk, where Ak is the stability matrix found from 
perturbations about the homogeneous state and which 
has entries given by Eq. l[2"3"]) . 

We are particularly interested in the situation where 
there is resonant behavior, that is, when there exist par- 
ticular frequencies when the denominators of Eqs. ([55]) 
and ([39]) arc small. The denominator vanishes when 
(iuj) 2 + {ioj) trAk + detAk = 0, which never occurs at 
real values of lo, however it does occur for complex to 
with non-zero real part if (trAk) 2 < 4 detAk- This pole 
in the complex lo— plane indicates the existence of a res- 
onance, and is exactly the same condition that the sta- 
bility matrix Ak has complex eigenvalues. This conforms 
with our intuition that the approach to the homogeneous 
stationary state needs to be oscillatory for demographic 
stochasticity to be able to turn this into cyclic behavior. 
If the lo dependence of the spectra numerators is ignored, 
then it is simple to show that the spectra have a maxi- 
mum in lo if additionally (trAk) 2 < 2 detAk- Using the 
full numerator results in a condition which is only slightly 
more complicated [H, G3 • 



B. Numerical results 

We used the stochastic simulations of the model de- 
fined by Eqs. ©-([I]) using the Gillespie algorithm [27l| . 
already mentioned in Section IIII CI to determine the 
Fourier transform of the fluctuations £fc(w) and r\k{Lo). 
These were then compared to those from the power spec- 
tra ((36]) and ([37]) . Once again we restricted ourselves to 
one dimension, which enabled us to obtain quite com- 
prehensive results. In practice the Fourier transforms 
are calculated by employing a discrete Fourier transform, 
and in order to compare the amplitudes obtained numeri- 
cally with the analytical results, the numerically averaged 
spectra contain an extra factor \ (iS x 5 t ) / (Af x J^t)\ 2 , where 
the 5 are the spacing between consecutive points and the 
Af the number of sampled points in space and time. 

We begin by showing the results of changing the num- 
ber of sites, 57. In Fig. [5] the left-hand column shows 
results obtained by taking = 200 with all other pa- 
rameters taking on the same values as in Section IIII CI 



The right-hand column shows results with the same pa- 
rameters again, except that = 500. The results from 
simulations were obtained by averaging 100 realizations 
of the process, taking an initial configuration to be the 
stationary state in the entire interval, and only once the 
oscillatory regime had been established. Specifically sim- 
ulation times were in the interval t e [1000, 2000]. 

The first two figures ([2]Ja) and (b)) show the typi- 
cal temporal evolution of the total population fractions. 
Subsequently, the results of simulations (upper graphs of 
Figs.[Hc)-(f)) and the analytic expressions ([55]) and ([3"9"]) 
(lower graphs of Figs. [U(c)-(f)) are displayed. Mention 
should be made of the scales of these (and subsequent) 
figures. The k take on discrete values 27m where n is an 
integer, since the length of the interval being considered 
is unity. In order to compare to the analytic forms, k is 
measured in units of 1/a, and so effectively it is ak which 
is plotted. This takes on discrete values 27m/f2, but we 
are looking at sufficiently large values of Vl that the k 
values appear continuous. For the lo— axis, the charac- 
teristic time which sets the scale is <5 t . It should also be 
noted that the k axis in Fig. [^Je) has been reversed to 
show the peak from another perspective. From Fig. [2jd) 
and (f), we see that the predator and prey spectra do 
not seem to differ appreciably. This was also found in 
the non-spatial case ll|. However, as we shall see later, 
if the migration rates are significantly different then the 
two spectra will differ. Also the fact that ak,n ^ 0, but 
that the analogous quantity in the non-spatial case, an, 
does vanish, leads to additional differences between the 
predator and prey spectra in the spatial version. 

For both values of O studied, we observe that the ana- 
lytic expressions and those obtained from simulating the 
full stochastic process show good agreement, which in- 
dicates that the use of the first two orders in the van 
Kampen approximation are sufficient for our purposes. 
We see that there is a large peak at a non-zero value 
of lo and so resonant behavior still occurs in this spatial 
model, just as it did in the non-spatial case. However, the 
height of the peak reduces with k and eventually at some 
finite value of k the peak disappears altogether. There is 
always an additional peak at lo — 0; this is much smaller 
and is just visible in Figs. [He) and (f). We will discuss 
it again shortly, when a different choice of the migration 
rates makes it far more prominent. 

In Fig. [3] similar plots are shown for two different val- 
ues of the migration rates \x\ and /Z2, keeping all other 
parameters as before (except in one case where we take 
c?2 7^ 0) and taking £1 = 500. The value of di was changed 
so that the fixed-point values <fi* and "0* were different, 
which made some of the plots clearer. 

Finally, as shown in Fig. [2J we found that making one 
migration rate considerably bigger than the other led to 
significant differences. Although the peaks at non-zero 
to were still present, they looked rather different for the 
predator and for the prey spectra. Also noteworthy is the 
peak at zero frequency, which is now much larger than 
before in the case of the prey. The graph is cut-off at 
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FIG. 2: Temporal evolution of the total population fractions and power spectra obtained from averaging 150 independent 
realizations with Q — 200 (left column), and averaging 100 realizations with a system composed by Q — 500 sites (right 
column). The reaction rates are the same as those indicated in Figure [T] The upper graphs in panels (c)-(f) show the results 
of the simulations while the lower graphs the analytic predictions (|38p ~ (|39p . 
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FIG. 3: Temporal evolution of the total population fractions and spectra obtained from numerical simulations of the process 
(upper graphs) and from Eqs. (138[) - (|39l ) (lower graphs). The site capacity and the number of sites were N = 500 and Q = 500. 
The left-column panels were obtained employing the same local reaction rates as in the previous cases and fii — 0.5Q,, \ii = 0.7S1, 
whereas the right-column panels were obtained with \i\ = 0.8H, [12 = 0.9Q, and di = 0.050. The spectra in both cases were 
obtained by averaging 100 independent realizations. 
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k ~ 1 only because it becomes much more noisy at larger 
values of k and so rather difficult to interpret. A similar 
result is obtained if we swap the values of the migration 
rates, but now it will be the predator fluctuations which 
will exhibit the large amplification effect. 



V. CONCLUSION 

In the work that we have presented here we have 
stressed the systematic nature of the procedures em- 
ployed and the generic nature of the results obtained. 
The starting point was the ILM (HJ)— (J3J) , but many of the 
results that we give are not sensitive to the precise form 
of the model employed. For instance, births and preda- 
tor events could have an alternative (or additional) rule 
which would involve nearest-neighbor patches. An exam- 
ple would be BiEj — > BiBj, where i and j are nearest 
neighbor sites, which would mean that a birth could only 
take place of there was space in the adjoining patch. The 
definition of the neighborhood could also vary to include 
next-nearest neighbors or a Moore neighborhood, rather 
than a von Neumann one. All these changes would give 
the same behavior at the population level, and in many 
cases exactly the same model, and leave the form of our 
results unchanged. 

In a similar way, the nature of the lattice, and its di- 
mension, only enter the differential equations through the 
discrete Laplacian operator Ak and factors of a d , leav- 
ing the essential aspects of quantities such as the power 
spectra unchanged. One consequence of this observation 
is that the very good agreement between the analyti- 
cally calculated power spectra and those found from the 
one-dimensional simulations should still occur in higher 
dimensions and for other models. This is the main justi- 
fication for restricting our simulations to one dimension 
and hence being able to obtain higher quality data. All 
these observations lead us to expect our results to be 
generally applicable and to be capable of straightforward 
generalization to other, similar, problems. 

The procedure we have followed is also systematic. 
Rather than writing down a PLM on phenomenological 
grounds, we have derived it within a expansion proce- 
dure with a small parameter (l/y/~N) from a more basic 
ILM. This allows us to relate the parameters of the PLM 
to those of the ILM, but also to derive the strength and 
nature of the noise that is a manifestation of the demo- 
graphic stochasticity, rather than putting it in by hand. 
The two sets of equations derived from the ILM — the 
macroscopic, or mean-field equations and the Langevin 
equations describing the stochastic fluctuations about the 
mean fields — capture the essential aspects of the dy- 
namics at the population level. Provided that N is not 
too small that stochastic extinction events are significant, 
they give a very good description of generic phenomena 
which one would expect to see in simple descriptions of 
systems with one predator species and one prey species. 

The main focus of this paper was on the power spectra. 



We found that the resonant amplification present in the 
well-mixed system is still present in the spatial system, 
although the height of the peak decreases with k, at least 
in the one-dimensional model. The spectra for the preda- 
tor and prey species can be made significantly different 
by making one of the migration rates much bigger than 
the other, a freedom that was not available to us in the 
non-spatial case. There is also a peak at u> = 0. This is 
present in the non-spatial model, but has no physical sig- 
nificance. Here it does: it corresponds to periodic spatial 
structures. This peak is very small if the migration rates 
are of the same order, but can be as large as the peak at 
lu =/= if the migration rates are sufficiently different. 

The existence of a large peak at non-zero u> and |k| 
means that when the system is studied at a spatial res- 
olution defined by k, there will large amplitude oscilla- 
tions of frequency cjo(k), where this is the position of 
the peak. While we can deduce the existence of such 
structures for general d from our analytic calculations, 
our numerical work has only been undertaken for d = 1. 
Since the topology of one-dimensional lattices constrain 
the dynamics from exhibiting more interesting structures 
in space and time (as have been reported in numerical 
studies of models of a similar nature 

EE EE El), these 

periodic structures may have more complicated forms in 
higher dimensions. 

The approach which consists of defining the time- 
evolution of a model by a master equation, and then 
performing some type of analysis which allows one to 
obtain not only the mean field theory, but corrections to 
it, has proved to be very effective in understan ding the 
results obtained from numerical simulations [HI, Il6l. |29|. 
In the case of the technique employed in this paper, there 
are many applications which can be envisaged — those 
which apply to completely different systems, but also 
predator-prey systems with a more complicated func- 
tional response. It would also be interesting to inves- 
tigate systems whose deterministic limit exhibits Turing 
instabilities [12, EH . In other words, the general approach 
we have discussed here, and the results we have reported, 
have a very general nature. This implies that resonant 
amplification of stochastic fluctuations will be frequently 
seen in lattice models and lead to cyclic behavior in a 
wide range of systems. 
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APPENDIX A: SYSTEM SIZE EXPANSION 

In this Appendix the master equation for the model 
discussed in the main text is expanded to leading order 
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FIG. 4: (a) Total population fractions and (b) spatial configurations for the predator and prey fractions, (c)-(d) Numerically and 
analytically obtained power spectra obtained from 70 realizations of the process and from expressions (|38|l and (|39[) respectively. 
The migration rates were /ii = 1.0Q, fi2 = O.Olfi, and the local rates are the same as in Figure [1] The amplification effect is 
stronger that in the previous cases particularly in the case of the prey spectra. Simulations have been carried out swapping the 
values of the rates, showing a similar effect, but for the other spectrum. 



(which gives the macroscopic laws) and next-to-leading 
order (which gives the linear noise approximations) in the 
van Kampen system-size expansion [6| . The system-size 
expansion is not usually applied to systems with spatial 
degrees of freedom (but see [30]), and there are a num- 
ber of possible ways of proceeding. Here we will take 
what is perhaps the simplest case, and assume that the 
expansion parameter is 1/yN, that is, each lattice site 
is treated as a subsystem for which the carrying capac- 
ity becomes large. The calculation may be performed in 
a way which is similar to the non-spatial case; whereas 
in the non-spatial model there were two degrees of free- 
dom: the number of predators, n, and the number of 
prey, m, there are now 2Q degrees of freedom, rn and m^, 
i = 1, . . . , f2. In what follows we will therefore limit our- 
selves to an outline of the method and to the statement 
of key intermediate results. For all fuller description of 
the method, reference should be made to van Kampen's 
book [6( or papers which apply the method to related 
problems [II [l3. 



The system-size expansion begins with the mapping 
^ = & + ^=lfc + (Al) 

Here <f>i(t) and ipi(t) will be the variables in the PLM, 
and the stochastic variables (t) and 77, (t) will appear in 
the Langevin equations at next to leading order. 

Under this transformation, the left-hand side of the 
master equation ([7]) becomes: 

where £j = —(N)irf>i, rji = —(N)^tpi and where n is 
the probability density function, but now expressed as a 
function of (j>i,if>i and t. To determine the form of the 
right-hand side of the master equation in terms of the 
new variables, we need to write 7J loc and 7^ mg , given by 
Eqs. © and (|T0|) respectively, in terms of these new vari- 
ables. This consists of two stages: first writing the step 
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operators @ as operators involving the new variables, 
and secondly, determining their action on the transition 
probabilities (|5|) and ([6]). 

Beginning with 7~ loc the first stage gives 



E r 



El 1 -! = 



E Xi E y, ~ 1 



dru 2 

, r 1 9 1 at 1 d2 
- N 5 «F + 53 + 



i^r- f- - 

2 % 



(A3) 



We can now list the various contributions we obtain, at 
order AT 2 and iV , which we need in order to find 7^ loc 
as defined in Eq. (J9j> : 

(i) (E Xi - 1) dinf. 



<9& 



N 



d_ 



(ii) (£^l)(^ + ( i 2 m,): 

1 d 
iV5 : (d 2 ipi + 2p2'ipi<t>i) 7r - 

d d 
N° : (d 2 + 2p 2 ^)o- Vi, 2p2Vi«~6. 



V'i +P2lpi(/>i 1 . , j 



( m ) W - 1) ( 



2bm,i(N —rii — mi) 

iv 



AT* : -26^(1-^-^) 

TV : 26(2^-l + 0i)^-7?i 
d 



d 2 

Zbipi—^, b%pi <j>i) —3 

(iv) (£-i^_l) 
AT* : -2pi^0i — , 2p 1 <j> i ip i — 

o& dm 

3d d 2 d 2 

N° : 2pic/>i—rii, 2p 1 tp l —^, pxcj)^—^, pifatpi-^, 
orji dru drif dQ 

d d d 2 

-2pi0i— r?i, :-2piipi— -2^0^ . 
9& <, 



Identifying the terms of order AT 2 on the right- and left- 
hand sides of the master equation gives the contributions 
of the local reactions to the macroscopic laws: 

-k = §&-7fV>i&, (A4) 

G?2 2p 2 2pi 

- -Vi (1 - Vi - ■ (A5) 

If a rescaled time, t = t/il, is introduced, then these 
equations are exactly the PLM of the non-spatial version 
of the model [ll|. This is as it should be, since without 
including the nearest-neighbor couplings in T™ 1S , the sys- 
tem is simply f2 copies of the non-spatial model. 

Performing a similar identification of both sides of the 
master equation, but now for terms of order N° gives a 
Fokker-Planck equation: 



^ = -EaoWC(*))n] + 5 Eac^ 



[Bairn , 

where we have introduced the notation £j = The 
function Ai(C) and the matrix Bij are given by 
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(A7) 



and 



Kjloc 
a ij. 11 



a ij,22 



nloc 



— (di<t>i + 2pii/>i<f>i) 5 t j , 

— (2bip l (1 - fa - ipi) + d 2 ip l 
2 (pi + P2) ipi<fri) Sij , 



1 



(A8) 



The superscript loc denotes their origin from the local 
reaction contribution of the master equation, and the 
subscripts 1 and 2 refer to fi = £ and (2 = t], respec- 
tively. These results agree with the non-spatial results 
found in [ll[, up to a factor of ft, as required. It should 
also be noted that the function Ai(C) is linear in £j and 
rji with coefficients which are exactly those which would 
be obtained from a linear stability analysis of Eqs. (|A4[) 
and (|A5j) j|. This is given in the main text by Eq. (|20|) . 
which agrees with the results in Eq. (|A7[) . By contrast 
the Bij cannot be obtained from the macroscopic results. 

Next we carry out the same procedures on the contri- 
bution due to migration, 7^ mlg . To do this, the operator 
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expressions listed below are required 
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dm 
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(A9) 



These operators possess the general structure N~zLi + 
N~ 1 L2 1 with L\ equal to a difference of first deriva- 
tives and 1/2 = L\/2. In addition the transition rates 
([51) have a common structure as functions of N which 

is p \ NFi + Fi + F3 + . . .V when written in terms of 

the new variables, with p = (zQ) or ^/(zQ), depend- 
ing on which term one is considering. The Fk depend on 
the macroscopic fractions (<pi and ipi) and on the stochas- 
tic variables (&), except for F% which only depends on 
the former. Therefore the form of the part of the master 
equation involving migration terms is 



N~'Li + 



N- L L 2 



P (nFi + 
+ F-lL 2 + 



N2F2 ■ 

n. 



n 

(A10) 



keeping only terms of the order required. This allows us 
to identify the three main contributions: 

(a) The order term is identified with the second 
term in the left-hand side of the master equation 
(Eq. (HH with ii = -(N)3(j)i and m = —(N)$ipi) 
which leads to 2Q independent macroscopic equa- 
tions. 

(b) The order N° term pL\F2 is of the same order as 
the time-derivative in Eq. (|A2|) . Since it involves 
only first-order derivatives in it will give contri- 
butions which will add to the Ai in Eq. (|A6|) found 
for the purely local terms in the master equation. 

(c) The order term pF\Li is also of the same order 
as the time-derivative in Eq. (|A2j) . Since it involves 
only second-order derivatives in £j it will give con- 
tributions which will add to the 0™ in Eq. (|A6|) 
found for the purely local terms in the master equa- 
tion. 



As an example, the term T ni +\ %n \\ nun in Eq. (J6]) 
when written out in the new variables gives 



Mi 



[{fa (1 - & - Vi)} N + {(1 — <fii — 



n. 



- fa (6 + Vi)} - & (& + m) 

In the notation we have introduced above 
Fi = <pj (1 - <fo - ipi) . 



(All) 



(A12) 



The second term in Eq. ([6]), T n ._ l/nj+1 \ n . n ., can be 
obtained from the first term by interchanging i and j 
(and this is still true when the operators are included in 
Eq. (JTUJ) ) , so adding these expression together we find 



2m 



E (fa -fa)+^2 (fa^i ~ fa^ 



(A13) 



To obtain this we have identified dH/d^i, for each i, with 
the corresponding term on the left-hand side of the mas- 
ter equation (|A2|) . Using the discrete Laplacian operator 



A/ i =-x;(/ J --/<) > 



(A14) 



this may be written as 



Mi 



(A15) 

A similar analysis may be carried out for the terms 



and 



( E yl E Vi 



\ E Vi E vl 



IT, 



mi + l.mj — l\mi,rrij , 



it. 



jni — l.rrij +1 \mi ,mj ■ 



This will give the same form as above, but with the obvi- 
ous changes /ii — > P2, ipi <-> <f>%, etc.. For the macroscopic 
contribution one thus finds 



M2 



[A^Pi + iPiAfa - (piAiPij 



(A16) 



Adding Eq. (|A15|1 to the right-hand side of Eq. (|A~4|) and 
Eq. (|A16|) to the right-hand side of Eq. (|A5|) gives the 
set of macroscopic laws Eqs. (fTTj) - (fT2"j) for each patch i. 

Returning to the stochastic contributions, the one of 
type (b) coming from the term 

( E Xi E Xj ~ l) E ni + i^ nj -i\ nunj , 



is the -FVtype term in Eq. (|A11|) . Explicitly this is equal 
to 



EL 
zfl 



E 



_d d_ 



[(l-0 i -^)^'-0ite + '7i)] n - 

(A17) 
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The term 



(e^e- 1 - 1) iVi, 



gives precisely the same contribution, and adding these 
together one finds 

[{A-^A + (A^)}6 



+ {(f>iA - (Afa)}rji 
This may be written as 



n. 



n 



(A18) 



(A19) 



M2 ^£r[A,2iCi + A,22J7i]n, 



(A21) 



where 

A,n = A-ViA+(A^i) , A,i2 = faA-(Afa) . (A20) 

In an analogous way, the migrational contributions from 
the third and fourth terms in Eq. ([6]) give (letting [i\ — > 

l"2, fa <-> V>i and & <-► r/j) 

A 

£1 ^ % 

where 

A,22 = A-faA+(Afa) , A,2i = ipiA-(Aipi) . (A22) 

The results (|A19[) - (| A22[) can also be obtained through 
a linear-stability analysis of the non-local terms in 
Eqs. (fTTj) - (fl"2")) . They represent diffusion and should be 
added to the terms in Eq. (|A7|) which represent reactions, 
to give the complete contribution in the first term on the 
right-hand side of the Fokker-Planck equation (|A6[) . 

Finally, there are the terms of type (c), which have the 
form piqZ/2. We have already discussed the F\ terms, 
and the operators may be read off from Eq. (j A9|) . 
The four terms corresponding to those in Eq. ^ are: 

r . _ 1 2 
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dr/j 
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(A23) 



In this paper we will only be interested in studying 
the equations satisfied by the stochastic variables Ci = 
(£») i = 1, • . . , fl, when the transients in the macro- 
scopic equations (fTTj) - (fT2")) have died away. Then fa and 
ipi are equal to their fixed point values fa and tp* respec- 
tively, which are independent of the site label i. Adding 
the four contributions (|A23p in this case gives 



z8i 



i,3 



9Q 



d 2 



n 



n. 



These contributions are diagonal in the predator-prey 
variables (there are no mixed derivatives involving £ 
and 77), but is not diagonal in the site variables (there 
are mixed derivatives involving i and j). Comparing 
Eq. (|A24[) with the Fokker-Planck equation (|A6|) . we see 
that the contributions to the matrix B, which add to 
those in Eq. (|A8|) are 



? mig _ 



R mig 



^•(1-0- 








(1 - fa 






-ip*)j(ij) 



(A25) 

where Juj) is zero unless i and j are nearest neighbors. 

In summary, the order ./V terms give the Fokker- 
Planck equation (|A6I) . with the function Ai(C) and the 
matrix Bij being given by: 



Ai,2 = Q!i,21^ + «i,22?7i , 



(A26) 



(A24) 



where the a are exactly the coefficients found in Section 
IIII Al by linear stability analysis, and 

B«,n = [{d l fa + 2p 1 ^*fa)+Aii l fa{l-fa-ip*))8 ij 

- —fa^-fa-tp*)^, 

% 22 = [(26^(1-^-^*) + ^* 

+ 2 (p! + pa) + W (1 - 0* - V*)] % 

- (1 - 0* - r) J m , 

By,12 = %21 = [-2pi0V1 % ■ (A27) 

In the above we have assumed that the Fokker-Planck 
equation (|A6j) has been re-expressed in terms of the 
rescaled time t = t/Sl, in order to eliminate factors of 
f2 _1 from A and B. 



APPENDIX B: FOURIER ANALYSIS 

As discussed in the main text we carry out a temporal 
Fourier transform in order to calculate the power spec- 
tra associated with the fluctuations about the stationary 
state in order to identify temporal cycles, but we also 
wish to carry out spatial Fourier transforms. There are 
a number of reasons for doing this: (a) the translational 
invariance of the stationary state means that quantities 
of interest become diagonal in Fourier space, (b) because 
of this the continuum limit is easily taken, and (c) the 
power spectra are naturally generalized from the non- 
spatial case to depend on the wave-vector as well as on 
the frequency. 
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We largely follow the conventions of Chaitin and 
Lubensky [3l| in introducing the spatial Fourier trans- 
forms. That is, we define the Fourier transform, / k , of a 
function /j defined on a d— dimensional hypercubic lat- 
tice, with lattice spacing a, by 



Now 



.ik.aj 



/k , 



(Bl) 



where, for clarity, we have deviated from the usual no- 
tation of the main text and written the lattice site label 
j as a vector. Here k is restricted to the first Brillouin 
zone: —(tt/o) < fc~ < (tt/o), 7 = 1, . . . , d. We will also 
require the result [3l[ 



-ik.aj = n5 . ( 



(B2) 



Using the definition (|B1[) we may take the Fourier 
transform of the Langevin equation (|28| . This is straight- 
forward for the time derivative on the left-hand side 
and for the noise term A;. For the At term we use 
Eq. (|A26|) where the a are made up of the local con- 
stant terms (|20)) and those coming from diffusion (|A20[) 
and (|A22j) . At the fixed point where 4> and ip are homoge- 
neous these diffusion operators are site-independent and 
given by D lx = (1 - 4>*)A, D 12 = <f>*A, D 2 i = i/)*A and 
D22 = (l — <t>*)A. The Fourier transform of the Langevin 
equation thus takes the form (|30[) . with the a given by 
Eq. (|2"3")l . where A k is the Fourier transform of the dis- 
crete Laplacian operator A. From the definitions (|A14[) 
and (jBip this is easily shown to be given by Eq. (|24|) . 

To complete the description of the Langevin equation 
in k— space, we need rewrite the correlation function (|29[) . 
Taking the Fourier transform of both Aj(r) and Aj(r') 
yields 

<A k (r)A k , (/)) = a 2d £ e~^ ai e-^'^ By 5{r - r') . 

' (B3) 
However, By is given by Eq. (j A27[) and is only non-zero 

if i = j or if i and j are nearest-neighbors. That is, it has 

the form 



B» = b<® <5n + & (1) J ( 



ij) 



(B4) 



The translational invariance of By is quite clear: it can 
be completely specified by the difference d = j i: 



B d = 6 (0) <5 d , + & (1) <5|d|,i 



(B5) 



Inserting the expression for Bd in terms of its Fourier 
transform, B q , in Eq. (|B3[) . we have from Eqs. (|B1[) and 
(|B2|) that 

(A k (r)A k ,(r')> = a d n ^ B q <5 k , q <5 k ,,_ q 5(r - r') 

q 

= B k a d nS k+ v t0 6{T-T'). (B6) 



B k = a d J2 e_ ' k ad b a 



cos (fc 7 a) 

.7=1 



(B7) 



using Eq. ([B5| . In terms of A k defined by Eq. ([M]), this 
may be written as 



Bv = a d 



(B8) 



since for a hypercubic lattice the coordination number is 
z = 2d. Writing these out explicitly using Eqs. (|A27[) 
and (IB4[) gives Eq. ([52^) in the main text. 

Finally, we can ask what happens as we take the lattice 
spacing, a, to zero, but keeping fla d (the area, if d = 2) 
fixed. Using Eq. ([54]) and 



cos (/c 7 a) ~ 1 



(A; 7 a) 



0((fca) 4 ) 



(B9) 



we see that A k 



2 k 2 /d + 0(k 4 ). Since A k always 



appears along with the migration rates, the factor of a 2 / d 
can always be absorbed into these rates by defining new 
quantities 



1 2 

Mi = -a Hi, M2 
d 



1 



a 2 n 2 ■ 



(BIO) 



So for instance, in Eqs. (|23|) and ([32)1 the A k can be 
replaced by —k 2 and fix and /12 by /ii and /i2 respectively, 
as a becomes small (or equivalently Q becomes large) . In 
this limit O a d Jic+ic^o becomes (27r) rf <5(k + k') [31I, and 
therefore Eq. ([3"Tj) becomes 



(A k (r)A k - (r')) = B k (27r) d ,5(k + k') S(r - r') , (Bll) 

where i? k is given by Eq. (|32[) . but with the small a ap- 
proximation described above. 

To obtain the power spectrum we need to take the 
temporal Fourier transform of Eq. (jBll[) . This yields 

(A k HA k , (c,/)) = B k (2Tr) d 5{k + k') (2tt) S (w + w') . 

(B12) 

Since there are only contributions in the above formula 
when k' = — k and u>' = -co this is frequently written as 

(A k (u;)A_ k (-u;)> = S k , (B13) 

or equivalently since A k (w) = A_ k (— u>), as in Eq. (J35J). 
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